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ABSTRACT 

We investigate the effects of gas-disk gravity on the planetesimal dynamics in inclined binary systems, where 
the circumprimary disk plane is tilted by a significant angle (ib) with respect to the binary disk plane. Our focus 
is on the Lidov-Kozai mechanism and the evolution of planetesimal eccentricity and inclination. Using both 
analytical and numerical methods, we find that, on one hand, the disk gravity generally narrows down the 
Kozai-on region, i.e., the Lidov-Kozai effect can be suppressed in certain parts of (or even the whole of) the 
disk, depending on various parameters. In the Kozai-off region, planetesimals would move on orbits close to 
the mid-plane of gas-disk, with the relative angle (/ ) following a small amplitude periodical oscillation. On 
the other hand, when we include the effects of disk gravity, we find that the Lidov-Kozai effect can operate 
even at arbitrarily low inclinations (is), although lower ib leads to a smaller Kozai-on region. Furthermore, 
in the Kozai-on region, most planetesimals' eccentricities can be excited to extremely high values (~ 1), and 
such extreme high eccentricities usually accompany orbital flipping, i.e., planetesimal orbit flips back and 
forth between anterograde and retrograde. Once a planetesimal reaches very high orbital eccentricity, gas drag 
damping will shrink the planetesimal orbit, forming a "hot planetesimal" on a near circular orbit very close to 
the primary star. Such a mechanism, if replacing the planetesimals and gas drag damping with Jupiters and 
tidal damping respectively, may lead to frequent production of hot- Jupiters.. 
Subject headings: Celestial mechanics - planetary systems: formation 



1. INTRODUCTION 

As of today, over 60 exoplanets have been found in binary 
star systems, and current observations show that the multi- 
plicity rate of the detected exoplanet host stars is aro und 17% 
jMugrauer & Neuhaused 120091: lEggenbergeJ 1201 Oh . Planet 
formation in binary system systems presents numerous chal- 
lenges, as each stage of the planet formation process can 
be affected by the binary companion. A crucial stage that 
may be particularly sensitive to binary effects is the accu- 
mu lation of 1-10 km-s ized planetesimals (see the review 
by iHaphiphipoud (1201 Oh and the references therein). Be- 
cause of the perturbations from the binary companion, plan- 
etesimals will be excited to orbits with high relative veloci- 
ties^, preventing or even ceasing their growth ( Heppenheimer 
[19781: iWhitmire et alj|1998h . In the past decade,' with sev- 
eral discoveries of exoplane ts in close binary of separation ~ 
20 AU (Oueloz et al. 20 00l: iHatzes et al.ll2003l: IZucker et alj 
I2004t [Correia et aL,2008t iChauvin et al.ff201 Ih . the issue of 
planetesimal growth in binary systems becomes more chal- 
lenging and therefore attracts many researchers as well as 
many dynamical and collision al stu dies (Marzari & Sc holll 
2000; Moriwaki & Nakagawai (2004t iT hebault et al. 2004 
2006, 20 08. ,2009^ .Thebaul 1201 Paardekoope r et aU 
2008; Sch oll et all 120071: iPaa rdekooper & Leinhardt l20Tr 
Kiev & Nelson' '2008'; 'Beauge et al. 2010; Giu ppone et al 
2011; Xie & Zhou 2008, 2009; Xie et al. 2010a,b). 

Most of previous studies had considered only coplanar or 
near-coplanar cases, where the tilted angle between the bi- 
nary orbital plane and the circumprimary disk plane was close 
to zero, i.e., ib ~ 0. In fact, the coplanar case is reasonable 
only if it is applied to relatively close binary systems with 
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separ ation less than ~ 40-200 AU (lHalelll994 iJensen et al] 
120041) . beyond which the distribution of is is likely to be ran- 
dom and therefore the highly inchned case is more relevant. 
Planetesimal dynamics in highly inclined bina ry systems have 
only been investigated by Marzari et alj ([2009b), and most re - 
cently (at the time of writing this paper ) by Xie et aU (1201 Ih . 
Fragner et al. (2011), and Batv gin et all (1201 Ih 

Marzari et aL (.200 9b) found that, due to the perturbations 
of a inclined binary companion, planetesimals' nodal lines be- 
came progressively randomized, raising their relative veloci- 
ties to the degree that planetesimal growth by mutual colli- 
sion was significantly prevented. Nevertheless, the gaseous 
protoplanetary disk was ignored in their study, where plan- 
etesimals were only subject to the gravity of the binary stars. 
In reality, the gaseous disk can generally have crucial effects 
on planetesimal dynamics through two factors. One is the 
hydrod ynamic dr a g forc e, which has been investigated in de- 
tail by IXie et alj ([2011J). When gas drag is included, it is 
found that planetesimals from the outer regions (where condi- 
tions are hostile to planetesimal accretion) jump inward into 
an accretion-friendly region and pile-up there. This is re- 
ferred to as the planetesimal jumping-piling effect (PJP), and 
its general result, as shown in Xie et al. (2011), is to form a 
severely truncated and dense planetesimal disk around the pri- 
mary, providing conditions which are favorable for planetes- 
imal growth and potentially allow for the subsequent forma- 
tion of planets. Another crucial factor is the gravity of gaseous 
disk, which has been studied recently by Fragner et al. (201 1) 
with a hydro-dynamical model and by Batvgin e t al.l ( 1201 1[) 
with an analytical model. Generally, it is found that the grav- 
ity could pull the planetesimals back towards the middle plane 
of gas-disk. With proper conditions, such as a massive gas 
disk and/or a large binary distance, Lidov-Kozai effect could 
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be suppressed regardless of ig- However. lFragner et alj ( 1201 Ih 
could only focus on several typical cases with a few plan- 
etesimals in relative short simulation timescale because of the 
large computational hours, while Batygin et al. (2011) only 
concentrated on cases of very wide binaries with separation 
of ~ 1000 AU, aiming to just identify the important physical 
processes at play. 

In this paper, we investigate the effects of gas-disk gravity 
on planetesimal dynamics in inclined binary systems through 
both analytical and numerical fashions. Analytically, we de- 
rived the condition at which Lidov-Kozai effect is turned off 
by the disk gravity. Numerically, we confirm our analytical 
results and provide a global quantitative view of the role of 
gas-disk gravity in a large parameter space. Furthermore, one 
specific attention is given to the role o f disk gravity in shap- 
ing the PJP effect found by iXie et alj {2011). The paper is 
outlined as follows. 

In section 2, we describe our disk model and the initial set 
up. In section 3, we analyze the secular motion of a planetesi- 
mal under the gravity of disk and stars, focusing on the evolu- 
tion of planetesimal inclination and eccentricity. The analyti- 
cal study is followed by the numerical simulations presented 
in section 4. In section 5, we discuss some issues, including 
the PJP effect, implication for hot-Jupiters and disk preces- 
sion. Then in Section 6 we present our summary. 

2. DISK MODEL 

In our model, planetesimals are assumed to be initially 
moving on circular orbits in the mid-plane0 of a gaseous disk 
around the central star with one solar mass (Mq). A compan- 
ion star with mass of Mb (a free parameter) is orbiting around 
the central star-disk system with an orbital semimajor axis of 
flfi (a free parameter), inclination of is (a free parameter, rel- 
ative to the mid-plane of the disk), eccentricity of eg = 
(constant). In this paper, for simplicity, we only consider the 
circular case (eg = 0) and focus on the effect of gas gravity. 
For the eccentric case eg + 0, the Lidov-Kozai mechanism it- 
selfis more complicated (Katz et al. 201 1; Lithwick & Naoi 
1201 ll) . and thus this case0is not addressed in this paper 

For the gas di s k, we use a 3-dimension steady model as in 
iTakeuchi & LinI (l2002h . In cylindrical coordinates (r, z), the 
disk density profile is 



Psir^z) ^ P()h\j^] exp 



2 \ 



Z 



and the gas rotation rate is 
£lg(r, z) = QK,inid 



2 / 2 \ 



(1) 



(2) 



where t^K.mid is the Keplerian rotation in the mid-plane, 
hg{r) - /io(r/AU)*'>'^-'-'^^ is the scale height of gas disk, fg is 
a scaling n umber wi t h resp ect to the minimum mass of so- 
lar nebulae dHavashil ( Il981h . MMSN hereafter), po = 2.83 x 
10"'"gcm"^ r = -0-5, ho = 0.33 x 10"^ and y6 is a free 

- This is equivalent to assume that their proper eccentricities (and inclina- 
tions) are equal to the forced one. 

If ffi = 0, the average Hamiltonian is axisymmetric, thus the verti- 
cal angular momentum is an integral of motion, and the planetesimal orbit 
can be well described with the classic Kozai effect. Otherwise, if eg > 0. 
the vertical angular momentum is not a constant any more, and the clas- 
sic Kozai effect should be modified with the so called Eccentic Kozai eff'ect 
(Uthwick & Naoz 2011.1 



parameter The surface density of the disk has a power-low 
form ofI.g = ^fhipQfgho(rlA\]f, where ;t = ^ + 1.25. Nom- 
inally, in this paper, we set A: = -1 as the stander case. The 
inner and outer boundaries of the disk are set as r^a =0.1 AU 
and Tout - 12.5 AU. Their values have little eff'ect on the final 
results as long as they are not very close to the planetesimals. 

Our disk model is a very simple one, which ignores the re- 
action of gas disk to the binary perturbations. In more real- 
istic situati ons, as shown in the simulations of Larwood et al] 
(i99_6^; IFragner & NelsonI ((2010). the disk will become ec- 
centric, develop a warp and precess under the perturbations 
of the comp a nion star. Nevertheless, as pointed out by 
IFragner et al.l (1201 lb (also see our discussion in section 5.2) 
, planetesimals' secular dynamical behaviors are similar both 
in the evolving and non-evolving disk models, and thus our 
choice of a steady model can be reasonable to at least a ze- 
roth order approximation. Furthermore, using such a simple 
gas disk model is much less time consuming in computing the 
disk gravity as compared to using a hydrodynamical code, al- 
lowing us to see the effect of disk gravity on a much longer 
timescale. In addition, our simple gas disk model is conve- 
nient for making some analytical studies. 

It also worthy noting that the gaseous disk would slowly 
relax to th e binary orbital plane on the viscous evolution 
timescale ( Fragner & Nelsonl2010l) . Thus the assumption of a 
constant and relatively large is in our model is only relevant if 
the viscous timescale, f,,,-, ~ l{ah^Q.\^ ^ii), is larger than the 
secular perturbation timescale, tsec ~ 2jt/B. Equating these 
two timescales, the critical viscous parameter can be derived 
as 



Therefore, a high inclined case, which studied in this paper, 
is relevant only for a < Uc- If otherwise, a > ac, it is likely 
to reduce to a near co planar case, which has been studied in 
many previous works (iMarzari & Schollll2000 |:lThebau lt et alJ 
l200aiPaardekooper et al.ll2008l:IXie & ZhouE oOS. 20^. 

3. ANALYTIC STUDY 

In this section, we analytically study the secular dynamics 
of a planetesimal under the gravitational perturbations from 
both the companion star and the disk. Our interests focus on 
the evolution of the planetesimal's orbital eccentricity and in- 
clination, aiming to see how the Lidov-Kozai effect operates 
if the disk gravity is included. For the sake of this derivation, 
we introduce two coordinate systems: (1) the disk coordinate, 
where the X Y -plane is set as the disk mid-plane with the X 
direction towards the ascending node of the binary orbit, and 
(2) the binary coordinate, where the XF-plane is set as the 
orbital plane of the binary star with the X direction the same 
as X . In the disk coordinate system, angular elements are 
marked with a superscript (" ' "). For example, / and Q de- 
note the orbital inclination and longitude of ascending node 
in the disk coordinate system respectively, while / and Q are 
those in the binary coordinate system. 

3.1. The Disturbing Function 

The disturbing function of the star-disk-planet system can 
be expressed as 

R^Rd+ Rb, (4) 

where Rd and Rg aie contributions from the gravity of disk 
and binary stars, respectively. 
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According to iNagasawa et alj (120001) (see the appendix of 
their paper), taking the second order approximation, Ro can 
be expressed as 



(5) 



where «, a, e and /' are the orbital mean frequency, semi- 
major axis, eccentricity, inclination (in disk coordinate) of the 
planetesimal. T and S are two characteristic frequencies (see 
the appendix of this paper for details of their definition and 
calculation) which, under the disk model assumed in section 
2, can be approximately fit by the following formulas. 



\k+\ 



r(fl) = /,x4.5xl0-'*(^) rad/yr. 



(6) 



S{a)^fgX 1.7 X 10 



Uu/ 



t+l/4 



rad/yr. 



Note, T is actually the apsidal recession rate of a planetesi- 
mal if the planetesimal is affected only by the disk gravity in 
the coplanar case (/ - 0). Our calculation of T is generally 
consistent with that of iBatvgin et al.l (1201 ll) (here = 1 cor- 
responds to a disk mass of ~ O.O2M0 in the figure 2 of their 
paper) who used a similar disk model but different comput- 
ing technics. However, we emphasize that T should be scaled 
with the local surface density as in equation ^ rather than 
with the total mass of the disk (as wa s done in figure 2 of 
iBatygin et al l (!2011l) and Ean. (30 ) in .Fragneret al.rd2011h ) 
Following Innane n et al.l (Il997h . the binary part of the dis- 



turbing function can be expressed as: 



Rg = ^-B[e^ - (I + 4e^ - 5e^ cos^ co) sin^ ij 



(7) 



where / and oj denote the orbital inclination and pericenter 
(in binary coordinate) of the planetesimal. The characteristic 
frequency B is actually the precession rate of the planetesimal 
caused by the secular binary perturbation in the coplanar case 
(i = 0), and in the first order it can be expressed as 



B ~ B, 



3GMb 



4«4(1- 4)3/2 



(8) 



However, such a first order expression can be rather inaccurate 
un less one u ses the second order correction (B2) as suggested 
bv lThebault et al., (,2006.) and .Giuppone et al. (.2011,) . 



B-Bi^Bi 



1 + 



32Mfi 



Mil 



4)^ 



(9) 



Hereafter, we adopt B - Bo if there is no specific explanation. 

We plot T, S , and B in Figure[T]for the standard case, where 
the companion has mass of Mb - O.SMq, semimajor axis of 
as = 50 AU, and the disk surface density slope ofk = - 1 . The 
blue doted line, red dashed line, and black solid line indicate 
T, S and Z? as a function of the semimajor axis of the plan- 
etesimal (a), respectively. As can be seen, S is much greater 
than T and B in the whole of the plotted region of the disk, 
while T is greater (less) than B in the inner (outer) region. We 
will show in the following subsections that such a picture of 
T, B and S determines the dynamical evolution of the plan- 
etesimal's orbit. 

3.2. Evolution of the Planetesimal Inclination 
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Fig. 1. — The values of T{a), S{a), and B{a). The blue dotted line and red 
dashed line indicate T(a) and S(a) of an MMSN disk with inner edge 0.1 AU 
and outer edge 12.5AU and surface density slope oi k = -1. The black solid 
line shows B(a) of a companion star at ag = 50 AU with O.SM©. 



As the planetesimal is initially moving on a circular orbit in 
the mid-plane of the disk, the initial e and / are approximately 
zero, thus we ignore quantities that are on an order of higher 
than o{e^), o(i 2) or o(ei ). The disturbing function relating to 
the inclination then can be reduced to 



R 



rd 



na 



\-Si'^ 



B sin^ 



(10) 



Considering the relation between /, / and and introduc- 
ing two new variables p - i' sin Q' and q - i' cos Q', then 
Lagrange's planetary equations (relating to / and Q ) can be 
written as (see the appendix for detailed derivation) 



dp B 

— --(B cos Hb + S)q -i — sin Hb, 

dt 2 

dq 
It 



■j-^[Bcos^iB + s)p, 



(11) 



where Ib is the angle between the disk plane and the binary 
orbital plane. Note S > B > and the initial condition po = 
qo - 0, thus the solution of p and q can be written as 



B sin(2/B) 



sin(/f), 



2/ 
B sm(2iB) 
IBcosUb + 2S 



[1 - COS(/f)] : 



(12) 



where / = -sj (B cos^ Ib + S)(B cos 2iB + S). The maximum 
value of /' (note that ; 



yjp^ + q^) is 
B sin(2;B) 
B cos 2/b + S 



(13) 



S and is as 



As S » B shown in Figure [T] thus / 
small as on an order of o(B/S). It means that the planetesimal 
will keep its orbital plane close to the disk mid-plane, having 
the relative titled angle / oscillating with a frequency of / ~ 
S and an amplitude of ~ B/S . Such an analytical result is 
consistent with the hydrodynamical simulation performed by 
lF?agner et aL (.201 1.) . which has shown that the disk gravity 
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would try to pull the planetesimal orbit back to the disk mid- 
plane, maintaining a small relative angle (see the figures 3 and 
10 in their paper). 

Recalling the approximation (quantities that are 0{e^), 
0(f) or 0(ei) or higher are ignored) adopted before our 
derivation, we thus emphasize that our analytical results about 
the evolution of planetesimal inclination remain valid only if 
the planetesimal eccentricity is not excited or remains at a low 
value. Such an assumption, however, will break down if the 
Lidov-Kozai eff'ect kicks in. In the following subsection, we 
will address this issue, deriving the conditions in which the 
Lidov-Kozai effect takes over and planetesimal eccentricity is 
excited. 

3.3. Evolution of the Planetesimal Eccentricity 

Following llnnanen et al.l (119971) . the Lagrange planetary 
equations describing the evolution of the planetesimal's or- 
bital eccentricity (e) and pericenter (w) can be written as 

de SB J 

e sin(2tij) sin (14) 

dt 2 

— ~B(2-5sm-LosiifiB] + D. (15) 
dt ^ ' 



Compared to the equation (5) in the paper of llnnanen et al.l 
([l997), here we add the term of contribution from the disk 
(D), ignore quantities that are on an order of o(e^) or higher 
because of the initial circular planetesimal orbit, and take / ~ 
Ib because /' is very small before e is excited according to 
equation (fTsT l. The disk contribution term (D) can be written 
as (see the appendix for the detail of derivation) 



D~ -T - 



S B cos^ is 
B cos Hb + S 



(16) 



Not43, as J = is a singular point in the Lagrange planetary 
equation, thus equation ( fTSI ) and ( fTSI ) cannot be applied to the 
case of Ib - 0. 

For the Lidov-Kozai effect to kick in, we expect da>/dt x Q. 
Using this condition to eliminate the variable co in equations 
( fTSl ) and ( fT6] l, we then have. 



de 
dt 



'5eB- 



'2B + D 
, 5B 



2B + D 



sin^iB — 



5B 



(17) 



In order to increase e, we need de/dt > 0, which then leads to 



< 2 + DIB < Ssin^/fi. 



(18) 



This is the condition for Lidov-Kozai eff'ect to operate under 
the gravity from both binary stars and the disk. For disk-free 
case, i.e., ZJ = 0, then equation (fTSl l is reduced to the classical 
one, i.e., Ib > arcsin( V2/5) ~ 39.2°. 

As D and B are functions of the semimajor axis (a), in- 
equation (fTST l actually produces two critical semimajor axes, 
a lower limit of a^i and an upper limit of flc2, which can be 
derived from 2 -i- D/B = and 2 -H D/B = 5 sin^ Ib, re- 
spectively. If the disk is not very tenuous, such as /g > 0.1, 
then 5 » B holds and thus equation ( fT6] l can be reduced to 
£) X - Bcos^ Ib- In such a case, we can solve Od and ad 

Setting i's = in equation.s (15) and (16), the binary and disk's contri- 
butions to ^ are 2B and -B - T respectively, which are obviously wrong, 
though their sum (B - T) is correct. 



analytically if A: - -1,Q 
aci~4.17AU 

ac2~4.17AU 



^ Mb. . 2. 
(sm iB + 1) 



-2/3 



1^(1 -4sin^.) 



-2/3 



,if Ib > 30° 



(20) 



Comparing to the classical disk-free case, where Lidov- 
Kozai effect takes place only if Ib > 39.2°, here Lidov- 
Kozai effect (or eccentricity excitation) can occur for an ar- 
bitrary Ib, and the value of Ib just determines the disk range 
(flci < a < ad) that subject to Lidov-Kozai effect. 

4. NUMERICAL STUDY 

In this section, we perform numerical simulations to test 
our analytical results presented in section 3. Planetesimals 
are only subject to the gravity from the binary stars and the 
dislS. We calculate the disk's gravity at lattice points in the 
r -z plane before orbital integrations and obtained the gravi- 
tational force at arbitrary point by bicubic interpolation (see 
the appendix for a detail description about computing disk 
gravity). The equations of planetesi mal motion are integr ated 
using a fourth-order Hermit method (iKokubo et al.lll998h . 

4.1. Examples 

As a first example (hereafter referred to as the standard 
case), we assume that Mb - O.SMq, as - 50 AU, Ib - 50° for 
the binary and fg = 1 , A: = - 1 for the disk. The results of this 
case is plotted in figure |2] and [3j 
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Fig. 2. — Planetesimal's maximum eccentricities (panel (a)) and inclinations 
(panel (b)) as a function of its initial semi-major axis in the stander case. Red 
dashed line indicates the results without including disk gravity. In panel (a), 
the vertical bla ck dash-dotted line indicates the analytical boundary of Kozai 
effect (Eqn ll9|. In panel (b), the black dashed line shows the analytical result 
from equation 1731 



In order to analytically derive a^i and with moderate accuracy, we 
make the compromise that B and B[ has the same dependency on a but a little 
difference in normalizati on, name ly B = (1 + r])B[. Here in equation )19t 

= 0.4. Note, equations )19t and 1201 cannot be applied to the disk free case 
by just setting fg = because we have presupposed that fg > 0.1. A nd f or 
the case of /: # -1, a^i and 0^2 should be solved numerically from Eqn IlSI 

* In fact, planetesimals are also subject to the hydrodynamical drag force 
from the gas disk. See section 5.1 for a discussion of gas drag or see the paper 
of Jfie et al. 1 2 01 1.) for a detailed study of the effects of gas drag. 
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Fig. 3. — Orbital evolution of a planetesimal with semimajor axis at 6.5 AU. All the orbital elements are in the binary coordinate. The two left panels are results 
of the standard case, while the two right panels are results for the case with the same binary configuration but without including disk gravity. (Note the different 
scales for the eccentricity and inchnation scales in the two plots.) 



1.0 
0.8 

0.6 
0.4 
0.2 
0.0 

100 

10 



0.1 



■ with disk 
without disk 
analytic u 




(a) 




4 6 8 10 

initial semi-major axis (AU) 



Fig. 4. — Similar to Figure|2]but with ig = 20° 



In Figure|2] we plot the maximum orbital eccentricity {e„ax) 
and inclination in disk coordinates) that the planetesi- 
mal achieved during its evolution as a function of its orbital 
semimajor axis. As can be seen from Figure |2] in the inner 
region, planetesimal eccentricities are not excited, and they 
remain at very low inclinations with fitting well with our 
analytical result (EqnfTSll. In the outer region, the Lidov- 
Kozai effect is switched on, and thus leads to large planetes- 
imal eccentricities (e,„aj: ~ 1) and inclinations (i,,,^^ > 90°). 
The boundary that separates the inner Kozai-off region and 
the outer Kozai-on region is roughly consistent with the ana- 
lytical estimate (EqnfT9]l. In addition, we also plot the results 
of the disk-free case as shown in the red dashed line in Figure 
|2] Comparing the two cases of with and without disk, we see 
that e,„av is much larger (close to unity) in the former case. 

In Figure|3] for a specific planetesimal with semimajor axis 
of 6.5 AU where the Lidov-Kozai effect should be switched 
on according to Figure |2] we plot the temporal evolution of its 
orbital eccentricity (e), inclination (i), longitude of periastron 
((L)) and ascending node (Q) for the two cases with and with- 
out disk. Note, here all the angular elements plotted in figure|3] 
are in the binary coordinate. In the case without disk, the two 



right panels show the classical "Lidov-Kozai" cycle where the 
eccentricity and inclination are evolving out of phase. How- 
ever, the situation is very different if the disk gravity is in- 
cluded in. In such a case, as shown in the two left panels of 
Figure|3] the planetesimal maintains its orbit around the initial 
one (/ = is, O - 180°) for a while at the beginning when the 
eccentricity is not very high. As the planetesimal eccentricity 
increases to the degree where e ~ I, the planetesimal quickly 
flips to a retrograde orbit but still in the same plane (the mid- 
plane of gas disk) with i ~ n - is and Q. - Qoin. We note 
that such an orbital flip as well as the associated high orbital 
eccentricity is ve ry simil ar to the one observe d recently by 
iNaoz etal.l (I2011 ab): Lit hwick & Naozl (12011b . where they 
assume a non-zero eccentricity of the outer perturbing body, 
and the orbital flip of the inner bo dy is due to the so-calle d 
Eccentric Lidov-Kozai Mechanism dLithwick & Naozll201 lb . 
While in the present paper, we assume a zero eccentricity of 
the outer perturbing body (eg - 0), and thus the orbital flip 
observed in Figure [3] should be due to the effect of disk grav- 
ity. 

As a second example (hereafter referred to as the low in- 
clination case), we just change the binary orbital inclination 
to is - 20° and keep all the other parameters the same as 
in the standard case. The result of this low inclination case 
is plotted in figure |4] In contrast to the is - 50° case, here 
the Lidov-Kozai effect can only take place within the region 
flci < a < flc2- This is consistent with our analytical results in 
equations [19] and |20] 

4.2. Parameter exploration 

In this subsection, we extend the standard case above by 
numerically investigating the effects of other parameters, in- 
cluding iB, as, Mb, fg, and k. We adopt the following strategy: 
To investigate the effect of a given parameter, we set this pa- 
rameter as the only free one and fix all other parameters as the 
same as in the standard case. The results are then plotted in 
Figure |5] and |6] The former shows the radial distribution of 
planetesimal's maximum eccentricity and its dependency on 
iB (top left), aB (bottom left), Mb (top right) and fg (bottom 
right) in the case of ^ = -1. The latter just shows the depen- 
dency on iB but for the cases of different k values. Some major 
features can be summarized as the following. 

As shown in Figure |5] (1) the Lidov-Kozai effect can be 
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switched on even with ig as small as ~ 5°, although the width 
of the Kozai- on region decreases as is decreases. (2) The 
Lidov-Kozai effect can be suppressed over a larger region if 
either the mass of the companion star decreases, or the separa- 
tion of the companion and/or the density of the disk increases. 
(3) The analytical results (dashed and solid lines, see also in 
EqnfT9landl20ll approximate the numerical results in the inner 
region with a < 9 - 10 AU. Beyond this, in the region close 
to the disk outer boundary and the orbital stability boundary, 
the deviation is large, indicating our analytical approximation 
is not valid there. And (4) the boundaries that separate the 
Kozai-on and Kozai-off regions are very steep; most planetes- 
imal eccentricities are either very high (close to 1) or very low 
(close to 0) with planetesimals of moderate eccentricities be- 
ing very rare. The effect of the disk density slope k can be seen 
from Figure |6] which shows (5) the Kozai-off region extends 
outwards more and more as k increases, i.e., the disk radial 
density profile becomes more flat. In the case of A; = -1/2, 
the Lidov-Kozai effect turns off in the whole disk. 

5. DISCUSSION 

5.1. Planetesimal Jumping and Pile-up (PJP) 

In the early stage, there must be a gas disk around the pri- 
mary star The gas disk has crucial effects on the dynam- 
ics of planetesimals through two factors. One is the gravity, 
which was studied in detail in previous sections of this paper 
The other one is the hydrodynam ic drag force, w hose role 
has been inves tigated in detail by Xie et al.l (1201 li) . In gen- 
eral, IxJe^eTari (j20 11) find that if planetesimals are excited to 
orbits with very high inclinations (relative to the disk plane) 
and eccentricities, they will be subjected to very strong hydro- 
dynamic drag forces from the gas disk, letting them jump in- 
ward and pile up, i.e., the so-called Planetesimal Jumping and 
Pile-up (PJP ) effect. Neverth eless, the disk gravity is not in- 
cluded in bv lXie et"aLl (1201 Ih. In the following, we show how 
the PJP effect is modified if both the gas drag and disk gravity 
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are included. 

We consider four cases, (a) the standard case as described in 
section 4.1, (b) a more compact case- similar to the standard 
case but with ag decreasing to 40 AU, (c) a low inclination 
case- similar to the standard case but with ig decreasing to 
20°, and (d) a disk gravity free case- similar to the standard 
case but the disk gravity is not included. In each case, gas 
drag force is calculated by assuming a single planetesimal ra- 



dial size of 5 km and following the procedure as described 
in section 2.2 ol 
figures|2]and[8] 



in section 2.2 of iXie et al.l (120 111) . The results are plotted in 



In Figure |7] we plot the evolution of orbital semimajor axis 
(also periastron, top panel) and eccentricity (bottom panel) of 
three planetesimals in case (a), i.e., the standard case. Three 
planetesimals with near circular orbits starting from 5, 6 and 
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Fig. 7. — Evolution of orbital semimajor axis (also periastron, top panel) and 
eccentricity (bottom panel) of three planetesimals with initial semimajor axes 
of 5 (green) , 6 (red) and 7 (blue) AU respectively. The binary configuration 
is set as in the standard case described in section 4. 1 but adding in the gas 
drag. The black dashed-dot line denotes the critical semimajor axis a^i = 5.4 
AU as measured from Figure |2] 
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Fig. 8. — Evolution of the local surface density enhancement (Z/Zq) of the 
planetesimal disk in the standard case (a), compact case (b), small inclination 
case (c) and the disk gravity free case (d). The vertical black dash-dotted 
lines show the critical semi-major axis that separate the Kozai-on and Kozai- 
off regions. The number on the left of each profile curve denotes the total 
number of planetesimals which have migrated to the innermost region within 
0.2 AU. Note the scale of the vertical axis is dilferent in the bottom panel as 
compared to those in the others. 

7 AU. Beyond = 5.4 AU, where the Lidov-Kozai effect 
is switched on, the two planetesimals' eccentricities are ex- 
cited and thus they suffer significant gas drag force, leading 
to rapid inward migration. The one starting from 7 AU is ex- 
cited to extremely high eccentricity and directly jumps to the 
innermost orbit, and the one starting from 6 AU with modest 
eccentricity quickly migrates to 2-3 AU the central star On 
the other side, for Oc < 5.4 AU, the Lidov-Kozai effect is sup- 
pressed, and thus the planetesimal starting from 5 AU does 
not suffer eccentricity excitation and hence does not migrate. 

In Figure |8] we plot the evolution of the local surface den- 
sity enhancement (2/Eo) of the planetesimal disk for the four 
cases. In each case, 5000 planetesimals treated as test parti- 
cle tracers, are initially distributed uniformly from 0.2 AU to 
10 AU in the mid-plane of gas disk with circular orbits, thus 
the initial profile follows a power law with a semi-major axis 
dependence equal to - 1 . By tracing the radial distribution of 
those planetesimals, we can calculate the local surface density 
enhancement (S/Sq) of the planetesimal disk. The results are 
shown in Figure [8] and can be summarized as follows: 

(1) As can be seen in the top panel (case a), planetesimals 
in the outer Kozai-on region migrate into the inner Kozai-off 
region (in fact, most are "jumping" as shown in Figure. |2]and 
pile up there, leading to surface density enhancement in the 
innermost region (a < 0.2 AU). This pile-up effect increases 
when the binary separation, gb, decreases (case b), because 
the outer Kozai-on region is larger and thus more planetesi- 



mals can move in and pile up. Conversely, when we reduce 
(case c), the pile-up effect is reduced, and the pile-up region 
shifts outward to 1- 2 AU, because the outer Kozai-on region 
shrinks and less particles are excited to eccentricity close to 1 
(see Figure. |4|i. 

(2) If the disk gravity is not included (c ase d), then the si tua- 
tion reduces to the situation considered in lXie et al.l (1201 Ih . In 
such a case, the Lidov-Kozai effect can only be suppressed by 
the gas damping and this only in the very inner region within 
1-2 AU, where gas density is high. Beyond 1-2 AU, planetes- 
imals experience the Lidov-Kozai effect and most (if not all) 
will migrate inward and pile up within ~ 0.2 - 1.0 AU, lead- 
ing to an average local den sity enhancemen t of S/Sq ~ 10 
(see also in the Figure. 9 of IXie et aTl ( 1201 Ih ). However, we 
note that there are many fewer planetesimals piling up within 
region < 0.2 AU in case (d) than in case (a). The reason is 
that the planetesimal eccentricity (see Figure. [3]) in case (d) 
is not high enough to let planetesimals directly jump into the 
innermost region < 0.2 AU. 

In a word, the role of the disk gravity playing in the PJP 
effect can be summarized as the following. On one hand, disk 
gravity reduces the average PJP effect because it reduces the 
Kozai-on region in the outer disk. However, on the other hand, 
the disk gravity significantly enhances the PJP effect in the 
innermost region (< 0.2 AU) as it increases the orbital ec- 
centricities of planetesimals in the Kozai-on region to values 
close to 1 . 

5.2. Effects of Planetesimal Collisions 

In this paper, the planetesimals are treated as test parti- 
cles and their mutual collisions are ignored. As planetesimals 
jump inward, their orbital eccentricities are very high and thus 
they are potentially subject to collisions of very high relative 
velocities, which can entirely disrupt themselves. To know 
how relevant the collisions could be, we estim a te the c olli- 
sional timescale (fcoi) first. Following IXie et aTl (1201 Obi) , fcoi 
in an inclined binary system can be estimated asQ. 

where /ice is solid density enhancement beyond the ice line, 
/?p is the planetesimal radii. Taking typical parameters, i.e., 
/g = 1-0,/jce = 4.2, Ib = 50°, Ma = Mq, /?p = 5 km, it gives 
fcoi ~ 8 X 10^(fl/AU)^ yr. As planetesimals complete their 
jumps typically in a timescale of lO'* - 10^ yr shown in Figure 
|7]and|8] thus we conclude that collisions have little effects be- 
fore or under the process of planetesimal jumping, but they do 
play important roles after plantesimals jumping inside 1 AU. 
Actually, t his ex pectation is confirmed by the simulations in 
IXie et al.) (1201 ll) . As shown in the figure 13 of their paper, 
in the first 10 yr, the collisional velocity is very high but the 
collisional frequency is rather low. Afterwards, collisions be- 
come more frequent as more planetesimals pile up in the inner 
region. At the same time, planetesimals are damped to near 
coplanar and circular orbits, leading to a friendly condition 
for subsequent planetesimal growth by mutual collisions. 

5.3. Implication to the Formation of Hot Jupiters 

The Lidov-Kozai effect induced by a companion star in a 
binary system has been suggested as an important mecha- 
nism for the formation of hot Jupiters dWu & Murray! 120031 : 

^ Combine Eqn (3), (6) and (7) in lXie et alj l2010bD 
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iFabrvckv & Tremainell200'7l) . If a planet's eccentricity is high 
enough that its periastron is very close to the star (say < 0. 1 
AU) during the Kozai cycle, then tidal dissipation can kick 
in, which may circularize and shrink the planet's orbit, finally 
letting it become a hot planet. However, to induce such a high 
eccentricity by the classical Lidov-Kozai eff'ect, it needs an 
extre mely misaligned confi guratiorfl (say > 85°, accord- 
ing to lWu & Murravl (120031) ). which is not common and thus 
lowers the chan ce of forming a hot- Jupiter. As estimated by 
IWu et al.l (l2007h . such a "stellar Kozai" mechanism can only 
produce 10% hot Jupiters. 

Nevertheless, the situation will be different if the disk grav- 
ity is included in the Lidov-Kozai effect. In such an case, 
almost in the whole Kozai-on region, the eccentricity can be 
excited to be an arbitrarily high value even with very low ini- 
tial binary inclination (is), which produces many more "hot 
planetsimals" (a < 0.2 AU) as shown in Figure |8] Similarly, 
if our model and results can be applied to a Jupiter-like planet 
(by replacing the gas drag damping with tidal damping in 
figure |8]l , it should also produce many more hot Jupiters. The 
key issue is, to what degree the production rate of hot jupiter 
can increase via the above "modified stellar Kozai" mecha- 
nism. We will address this in detail in a forthcoming paper. 

5.4. Disk Precession 

In this paper, we assume that the gas disk is non-evolving 
and axisymmetric, which is apparently a crude approxima- 
tion. In fact, the gas disk (if it is not entir ely disrupted) should 
undergo a near rigid bo dy precession CLarwood et al.l[T996t 
iFragner & Nelsonll2010l) . and the precession rate can be esti- 
mated as 



eccentricity 



3GMb 



■ cos Ig 



(22) 



For the standard case considered in this paper, equation |22] 
gives Qrf 2.6 x 10"'*radyr"^ Adding such a rigid preces- 
sion to the gas disk, we re-run the simulations shown in the 
top-left panel of Figure |5] and plot the results in Figure|9] The 
two black solid curves in Figure |9] are two critical semimajor 
axes (flci and 0^2) derived from equation [TS] (not from EqnfT9l 
andl20]i by assuming 



D kT + B cos is + cos ;bO, 
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Fig. 9. — Panel (a). Maximum eccentricities distribution in tlie plane of 
if the precession of tlie disk is considered. Tlie black solid lines indicate the 
analytic boundaries of Kozai-on region. Panel (b). Same as the top left panel 
of Fig 5. Maximum eccentricities distribution if the precession of the disk is 
not considered. 



The critical ig can be lower if one considers the effect of the binary 
eccentricity iLithwick & Naoz 2011) 

' Note, the results might be different because the giant planet can signifi- 
cantly affect the gas disk, e.g., opening a gap. 
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Fig. 10. — Comparison between our results (left three panels) and those in 
figure 10 of Fragner et al. 12011) (right three panels). From top to bottom, 
they are evolutions of planetesimals' eccentricities, orbital inclinations (rela- 
tive to the binary orbital plane), and nodal precession, respectively. The line 
color indicate the planetesimal's initial semi-major axis. The initial setup, 
incl uding the configurati on of the binary and disk is adopted from the model 
3 of lFragner et al.l 1201 ID (see section 3.2 and table 2 of their paper). 



Although equation |23] is a very crude approximation, it pro- 
duces reasonable a^i and a^o which fit the numerical results 
as well as shown in Figure |9] Furthermore, both the cases 
with and without disk precession (comparing the two pan- 
els of Fig|9j produce some similar features, such as: (i) in 
the central regions of the disk, the Lidov-Kozai effect can be 
switched on at very low inclinations, and (ii) once the Lidov- 
Kozai effect is switched on, the planetesimal eccentricities can 
be much higher (most are close to 1) than those in the case 
without disk gravity. 

5.5. Comparison to the Hydrodynamical Results 

In order to further examine the validity of our disk model, 
we compare th e results of our mode l to the hydrodynamical 
results given bv IFragner et al.l (1201 lb . We adopt the same ini- 
tial set up as in the simulations shown in the figure 10 of 
IFragner e t al.l ( 1201 1|) and run the simulation with our model 
and numerical method described in Appendix D. The compar- 
ison results are plotted in Figure [10] As can be seen, the re- 
sults computed by our mod el are generally con sistent with the 
hydrodynamical results of lFragner et al.l ( 1201 1). Given such 
a comparison, we then feel confident of the results shown in 
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other places of this paper. 

6. SUMMARY 

In this paper, we investigated the effects of gas-disk gravity 
on planetesimal dynamics in inclined binary systems using 
both analytic and numerical methods. Our major conclusions 
are summarized as the following. 

Analytically, we derive that the planetesimal inclination fol- 
lows a small amplitude oscillation around the mid-plane of 
disk (see Eqn[T2] and [TsT l if the Lidov-Kozai effect is sup- 
pressed and thus planetesimal eccentricity is not excited. Fur- 
thermore, we derive the threshold condition (see EqnlTS] [T9l 
andl20]i in which the Lidov-Kozai effect switches on. We find 
the Lidov-Kozai effect can operate at very low inclinations if 
the disk gravity is considered, although the radial extent of the 
Kozai-on region is much smaller. 

Numerically, we confirm our analytical results over a very 
large parameter space by considering the variation of ig, as, 
Mb, fg- We find that the disk gravity narrows down the Kozai- 
on region, but at the same time significantly increases the 
maximum eccentricity (close to 1) of planetesimals in the 



Kozai-on region (see Figure |2|. Such high planetesimal ec- 
centricities usually accompany orbital flipping (see Figure[3l), 
i.e., planetesimal orbits flip back and forth between prograde 
to retrograde. 

Applying the effects of disk gravity to the planetesimal 
jumping-piling (PJP) process. We find that, on the average 
over the disk, disk gravity reduces the PJP effect. However, 
PJP effect is significantly enhanced in the innermost region 
within 0.2 AU (see Figure[8]) . In addition, given the extremely 
high eccentricity under the effects of disk gravity, we believe 
that the production rate of hot-Jupiters via the "stellar Kozai" 
mechanism could be increased. 
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APPENDIX 

A. THE DISTURBING FUNCTION OF THE DISK 



According to Nagasawa et al. 2000, taken to second order in e and the disturbing function caused by the disk can be 
expressed as 

Rd = [T{a)e^ +S(a)i'^] , (Al) 



T(a) and S (a) are given using an integral of cylindrical coordinates (r', 4>',z'): 



3(a - r' cos0')^ 



A5 



A3 



Gpg{r',zydr'd(l>'dz', 



(A2) 



r- r r p 

2n Jr,„ J-oo Jo \ 
where A - (a^ + r'^ + z'^ - 2ar' cos ^'^^^ 



■cos^'/fl 3z'2\ , 
-^3 I CrPg(f ,z)r dr d(p dz , 



B. INCLINATION EVOLUTION EQUATION 

Ignoring and higher order terms in the disturbing function R, perturbation function relating to the incUnation has the form 

-,2 



7? = — [-S(a)i'^ - B{a) sin^ i\ . 



(Bl) 



where / is the inclination in the binary coordinate and is that in the disk coordinate. In the binary coordinate system, the xy- 
plane is the binarys orbital plane, and the x-axis is the ascending node of the companion with respect to the disk. In the coordinate 
system of the disk, the x-axis is same as that of the binary coordinate system, and the xy-plane is the mid-plane of the disk. 
According to the geometrical relationship between the two coordinate, we have 



sin/ sin Q ^ 
- sin / cos Q 
cos/ 



1 
cos is sin is 
- sin is cos is ) 



sin /' sin Q' 
- sin /' cos Q' 
cos /' 



(B2) 



It is easy to obtain 



sin / sin Q = sin /' sin Q' , 

sin / cos Q = sin /' cos Q! cos is - cos /' sin is, 

cos / = sin /' cos Q' sin is + cos /' cos is- 



(B3) 
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Then 

sin^ i = l - sin^ i' cos^ Q' sin^ - cos^ i' cos^ (b - 2 sin i' cos CI' sin ig cos i' cos j'^ 
= sin^ is + sin^ cos^ - sin^ cos^ Q.' sin^ - sin i' cos cos Q' sin2jB 
= sin^ iB + sin^ sin^ Q' cos^ is + sin^ /' cos^ Q' cos Hb - sin /' cos i' cos Q' sin lis 
- i'^ sin^ D.' cos^ /s + /'^ cos^ Q' cos 2/b - i' cos Q' sin 2jb + sin^ Ib + o(j'^). (B4) 
If we ignore f' and higher order terms, the relationship becomes 

sin / - p cos iB + q cos lis - q sin(2/B) + sin Ib, (B5) 
where p = i' sin D.' and q = i' cos O'. Thus the perturbation function becomes 

2 

/? = ^ [-(fi cos^ iB + S)p^ - (S + B cos 2js)g2 + B smi2iB)q - B sin^ Ib] (B6) 
Using Lagrange's equations of motion, the evolution of the inclination is given by 

^ = - (Z? cos 2iB + S)q + ^ sin(2i£), 
at 2 

^=(Bcos'^iB + s)p. (B7) 
at ^ ' 



C. THE CONTRIBUTION OF DISK TO THE PERIASTRON PRECESSION 
The disturbing function of the disk has the form 

2 

Using Lagrange's equations of motion, and ignoring the term, we can the expression for the evolution of (x> cause by the disk 



na r t n^ 
Rd^ — — [T{a)e^ + S (a)i'^\ . (CI) 



da)\ 



-T+-S cotm'^/di) = D. (C2) 
dt I disk 2 

where i and i' are the inclinations in the companion coordinate system and disk coordinate system. Proceeding with the same 
method as used in Appendix A, we have 

sin^ i' = sin^ i sin^ Q cos^ Ib + sin^ i cos^ Q cos 2iB + sin i cos i cos Q sin 2iB + sin^ Ib, (C3) 
then we can obtain that 

5(sin'^ i') I di = sin 2/ sin^ CI cos^ iB + sin 2; cos^ Q cos 2;^ + cos 2/ cos Q sin 2/^ 

sin / sin Q cos is + sin / cos Q cos 2/^ + sin ; cos ; cos CI sm 2iB + sm ii, \ — cos CI sin 2iB — 2 cot i sm ii, 
= [2 sin^ i' cos i - 2 sin is (sin ; cos CI cos j'b + cos j sin Ja)J / sin j 

= {2 sin^ j' cos J - 2 sin i' cos Q' sin j'^) / sin i. (C4) 
Because initially i' = 0,we ignore o{i'^) term and have 

5(j"^)/5j = -2i' cos Q' = -2q. (C5) 

We have obtained previously that 

q = i' cos CI' = [1 - cos(/0] . (C6) 

2ccos2jb + 25 

For the case S » B, the timescale of the evolution of i' and CI' is much shorter than the Kozai timescale. Thus we replace q with 
its average value 

Bsin(2;fl) 
2Bcos2iB + 2S 

and D becomes 

r. T ^c- f T SBcos^iB 

D = -T--Si„axCoti = -T-— ■ (C8) 

2 (Bcos2jb+5) 
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D. GRAVITATIONAL FORCE OF THE DISK 

According to Nagasawa et al. (2000), the potential of the disk at (r, (p, z) is 

-In 



V^G T" r r p(r',z'ydcf>'dz'dr' 

X,„ J-ooJo (r2 + r'2-2rr'cos^^' + (z-z')2 + e)i/2' 

where rin, and rout are the inner edge and the outer edge of the disk, respectively, and e is a softening parameter used to avoid 
a singularity. For reasons of efficiency and precision, we set it to be 1 x 10"'. Derivative of the potential with respect to r or z 
yields the r or z component of the disks gravity. 



2n 



p(r',z'){r - r' cos (p')r'd(p'dz'dr' 



X,„ j-ooJo (r2 + r'2-2rr'cos0' + (z-z'P + e)^^2' 
F = r r r P(r',z')(z-z'ydcf>'dz'dr' 

' j,„ j-ooJo (r2 + r'2 - 2rr' cos 0' + (z-z')2 + 6)3/2' ^ > 

We numerically integrated equation (ID2I) using closed Newton-Cotes formulas with Bodes rule jPress et alj|1992l) . Since the 
integration costs too much CPU time, we can not do it for each orbital integration step. Instead, we calculated the disks gravity at 
lattice points in the r-z plane before starting the orbital integrations and obtained the gravitational force at arbitrary points during 
the orbital integration by performing bicubic interpolations (Press et al. 1992) using the value at lattice points. 
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